#Importing libraries for Shiny dashboard and documentation
library(tidyverse)
## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr 1.1.4 ✔ readr 2.1.5
## ✔ forcats 1.0.0 ✔ stringr 1.5.1
## ✔ ggplot2 3.5.1 ✔ tibble 3.2.1
## ✔ lubridate 1.9.4 ✔ tidyr 1.3.1
## ✔ purrr 1.0.2
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag() masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors
library(rlang)
##
## Attaching package: 'rlang'
##
## The following objects are masked from 'package:purrr':
##
## %@%, flatten, flatten_chr, flatten_dbl, flatten_int, flatten_lgl,
## flatten_raw, invoke, splice
library(roxygen2)
library(gridExtra)
##
## Attaching package: 'gridExtra'
##
## The following object is masked from 'package:dplyr':
##
## combine
library(here)
## here() starts at /Users/kaylaemerson/Desktop/ENV710_Statistics/FinalProject/ENV710_FinalProject
library(moments)
library(dunn.test)
library(nlme)
##
## Attaching package: 'nlme'
##
## The following object is masked from 'package:dplyr':
##
## collapse
library(gt)
#Load dataset from an RDS file named "optdata.rds"
settlement_data <- read_csv("Settlement_Data_with_updated_demographics.csv")
## New names:
## Rows: 411 Columns: 16
## ── Column specification
## ──────────────────────────────────────────────────────── Delimiter: "," chr
## (4): MonitoringYear, SettlementName, MPAName, TimePoint dbl (12): ...1,
## SettlementID, MPAID, Treatment, MAIndex, MTIndex, FSIndex, P...
## ℹ Use `spec()` to retrieve the full column specification for this data. ℹ
## Specify the column types or set `show_col_types = FALSE` to quiet this message.
## • `` -> `...1`
head(settlement_data)
#Rename columns
settlement_data <- settlement_data %>%
rename(`Household_Material_Assets` = MAIndex, `Household_Marine_Tenure` = MTIndex,
`Food_Security_Index` = FSIndex, `Place_Attachment_Index` = PAIndex,
`School_Enrollment_Rate` = SERate)
head(settlement_data)
#Round values to 2 decimal place
settlement_data <- settlement_data %>%
mutate(`Household_Material_Assets` = round(`Household_Material_Assets`, 2)) %>%
mutate(`Household_Marine_Tenure` = round(`Household_Marine_Tenure`, 2)) %>%
mutate(`Food_Security_Index` = round(`Food_Security_Index`, 2)) %>%
mutate(`Place_Attachment_Index` = round(`Place_Attachment_Index`, 2)) %>%
mutate(`School_Enrollment_Rate` = round(`School_Enrollment_Rate`, 2)) %>%
mutate(`YrResident` = round(`YrResident`, 2)) %>%
mutate(`Age`= round(`Age`, 2)) %>%
mutate(`Gender` = round(`Gender`, 2))
#Look at summary
summary(settlement_data)
## ...1 SettlementID MonitoringYear SettlementName
## Min. : 1.0 Min. : 1.00 Length:411 Length:411
## 1st Qu.:103.5 1st Qu.: 28.00 Class :character Class :character
## Median :206.0 Median : 54.00 Mode :character Mode :character
## Mean :206.0 Mean : 54.65
## 3rd Qu.:308.5 3rd Qu.: 80.00
## Max. :411.0 Max. :115.00
## MPAID MPAName Treatment TimePoint
## Min. :1.000 Length:411 Min. :0.0000 Length:411
## 1st Qu.:2.000 Class :character 1st Qu.:0.0000 Class :character
## Median :3.000 Mode :character Median :1.0000 Mode :character
## Mean :3.416 Mean :0.6886
## 3rd Qu.:5.000 3rd Qu.:1.0000
## Max. :6.000 Max. :1.0000
## Household_Material_Assets Household_Marine_Tenure Food_Security_Index
## Min. :0.1500 Min. :1.000 Min. :0.420
## 1st Qu.:0.4300 1st Qu.:1.900 1st Qu.:3.075
## Median :0.4900 Median :2.340 Median :3.480
## Mean :0.4821 Mean :2.401 Mean :3.406
## 3rd Qu.:0.5400 3rd Qu.:2.875 3rd Qu.:3.820
## Max. :0.7400 Max. :4.250 Max. :4.850
## Place_Attachment_Index School_Enrollment_Rate YrResident Age
## Min. :3.430 Min. :0.1200 Min. : 1.53 Min. :28.40
## 1st Qu.:3.900 1st Qu.:0.7400 1st Qu.:24.24 1st Qu.:41.84
## Median :4.000 Median :0.8300 Median :30.30 Median :44.82
## Mean :4.015 Mean :0.8123 Mean :29.06 Mean :44.69
## 3rd Qu.:4.090 3rd Qu.:0.9100 3rd Qu.:35.38 3rd Qu.:47.65
## Max. :4.990 Max. :1.0000 Max. :46.67 Max. :59.80
## Gender
## Min. :0.6000
## 1st Qu.:0.8800
## Median :0.9200
## Mean :0.9144
## 3rd Qu.:1.0000
## Max. :1.0000
#Remove NAs
settlement_data <- na.omit(settlement_data)
# Graph 1
# visualize the distribution of the food security index observations by MPA
histogram_1 <- ggplot(settlement_data, aes(x = Place_Attachment_Index,
fill = MPAName)) +
geom_histogram(bins = 8) + # follow rice's rule to get 6ish
labs(x = "Place Attachment Index",
y = "Settlement Observation Count",
title = "Place Attachment Index Distributions within Six MPAs in Papua, Indonesia") +
facet_wrap(~MPAName, nrow = 2)+
theme_bw() +
theme(legend.position = "none",
plot.title = element_text(hjust = 0.5))
# view histogram
histogram_1
# save histogram to files
ggsave("BarfieldAEmersonKSwarupSFig1.jpg", plot = histogram_1)
## Saving 7 x 5 in image
# Graph 2
# A scatterplot depicting the relationship between the food security
# and place attachment indices
fs_pa_scatterplot <- ggplot(settlement_data, aes(x = YrResident,
y = Place_Attachment_Index)) +
geom_point(aes(colour = MPAName)) +
geom_smooth(method = "lm", color = "black") +
facet_grid(~MPAName) +
labs(x = "Household Years of Residency",
y = "Place Attachment Index",
color = "MPA Name") +
ggtitle("Relationship Between Household Years of Residency
and MPA Place Attachment within Papua, Indonesian") +
theme_bw() +
theme(legend.position = "none",
legend.background = element_rect(color = "black"),
strip.text.x = element_text(size = 7),
plot.title = element_text(hjust = 0.5))
# View scatterplot
fs_pa_scatterplot
## `geom_smooth()` using formula = 'y ~ x'
A series of scatterplots depicting the relationship between food security and place attachment within specific marine protected area (MPA) regions in Indonesia. Both indices are scaled 1-5.
# save scatterplot to files
ggsave("BarfieldAEmersonKSwarupSFig2.jpg", plot = fs_pa_scatterplot)
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
#Household Marine Tenure is defined as the average number of rights a household
# has over marine resources in their MPA
#Potential for this to be a strong predictor of Place Attachment Index
#Draw a scatterplot with a regression line to examine correlation
mt_pa_scatterplot <- ggplot(settlement_data, aes(x = Household_Marine_Tenure,
y = Place_Attachment_Index)) +
geom_point() +
geom_smooth(method = "lm", col = "red") +
labs(x = "Household Marine Tenure",
y = "Place Attachment Index") +
ggtitle("Relationship of Place Attachment to Bird's Head Seascape MPAs against
Household Rights over Marine Resources")
mt_pa_scatterplot
## `geom_smooth()` using formula = 'y ~ x'
A scatterplot depicting the relationship between household marine tenure and place attachment over all marine protected areas (MPAs) in Bird’s Head Seascape, Indonesia. Both indices are scaled 1-5.
# save scatterplot to files
ggsave("BarfieldAEmersonKSwarupSFig3.jpg", plot = mt_pa_scatterplot)
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
Figure 4 (Sameer)
#Homoscedastic distribution of Place Attachment Index values when plotted against
# Household Marine Tenure in Figure 3
#Will examine scatterplot of Place Attachment Index against Household Marine
#Tenure for the various MPAs
mt_pa_all_mpas_scatterplot <- ggplot(settlement_data, aes(x = Household_Marine_Tenure,
y = Place_Attachment_Index)) +
geom_point(aes(colour = MPAName)) +
geom_smooth(method = "lm", color = "black") +
facet_grid(~MPAName) +
labs(x = "Household Marine Tenure",
y = "Place Attachment Index",
color = "MPA Name") +
ggtitle("Relationship Between Household Marine Tenure and Place Attachment to
MPAs in Bird's Head Seascape, Indonesia") +
theme_bw() +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black"),
strip.text.x = element_text(size = 7),
plot.title = element_text(hjust = 0.5))
#View scatterplot
mt_pa_all_mpas_scatterplot
## `geom_smooth()` using formula = 'y ~ x'
A series of scatterplots depicting the relationship between household marine tenure and place attachment within specific marine protected area (MPA) regions in Indonesia. Both indices are scaled 1-5.
# save scatterplot to files
ggsave("BarfieldAEmersonKSwarupSFig4.jpg", plot = mt_pa_all_mpas_scatterplot)
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
#Household Material Assets is defined as the economic well-being of a household
#Potential for this to be a strong predictor of Place Attachment Index
#Draw a scatterplot with a regression line to examine correlation
ma_pa_scatterplot <- ggplot(settlement_data, aes(x = Household_Material_Assets, y = Place_Attachment_Index)) +
geom_point() +
geom_smooth(method = "lm", col = "red") +
labs(x = "Household Material Assets",
y = "Place Attachment Index") +
ggtitle("Relationship of Place Attachment to Bird's Head Seascape MPAs against
Household Material Assets")
ma_pa_scatterplot
## `geom_smooth()` using formula = 'y ~ x'
A scatterplot depicting the relationship between household material assets and place attachment over all marine protected areas (MPAs) in Bird’s Head Seascape, Indonesia.
# save scatterplot to files
ggsave("BarfieldAEmersonKSwarupSFig5.jpg", plot = ma_pa_scatterplot)
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
Figure 6 (Sameer)
#Homoscedastic distribution of Place Attachment Index values when plotted against
# Household Material Assets in Figure 5
#Will examine scatterplot of Place Attachment Index against Household Material Assets for the various MPAs
ma_pa_all_mpas_scatterplot <- ggplot(settlement_data, aes(x = Household_Material_Assets,
y = Place_Attachment_Index)) +
geom_point(aes(colour = MPAName)) +
geom_smooth(method = "lm", color = "black") +
facet_grid(~MPAName) +
labs(x = "Household Material Assets",
y = "Place Attachment Index",
color = "MPA Name") +
ggtitle("Relationship Between Household Material Assets and Place Attachment
to MPAs in Bird's Head Seascape, Indonesia") +
theme_bw() +
theme(legend.position = "bottom",
legend.background = element_rect(color = "black"),
strip.text.x = element_text(size = 7))
#View scatterplot
ma_pa_all_mpas_scatterplot
## `geom_smooth()` using formula = 'y ~ x'
A series of scatterplots depicting the relationship between household material assets and place attachment within specific marine protected area (MPA) regions in Indonesia.
# save scatterplot to files
ggsave("BarfieldAEmersonKSwarupSFig6.jpg", plot = ma_pa_all_mpas_scatterplot)
## Saving 7 x 5 in image
## `geom_smooth()` using formula = 'y ~ x'
# Define Research Question:
# Is there a relationship between place attachment index and gender,
# age, and years of residency in Indonesian households?
# Examine data and possible correlations
# Raw Counts
ggplot(data = settlement_data, aes(x = Place_Attachment_Index)) +
geom_bar() +
labs(y = "Count") +
theme_bw()
ggplot(data = settlement_data, aes(x = Age)) +
geom_histogram() +
labs(y = "Count") +
theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = settlement_data, aes(x = Gender)) +
geom_bar() +
labs(y = "Count") +
theme_bw()
ggplot(data = settlement_data, aes(x = YrResident)) +
geom_histogram() +
labs(y = "Count") +
theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Possible Correlations
ggplot(data = settlement_data, aes(x = Age, y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = Gender, y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = YrResident, y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
# Multi-collinearity tests
cor.test(settlement_data$Age, settlement_data$Gender)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Age and settlement_data$Gender
## t = -6.1839, df = 409, p-value = 1.517e-09
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3784339 -0.2013749
## sample estimates:
## cor
## -0.2924084
cor.test(settlement_data$Age, settlement_data$YrResident)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Age and settlement_data$YrResident
## t = 9.553, df = 409, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3446207 0.5030580
## sample estimates:
## cor
## 0.4271122
cor.test(settlement_data$Gender, settlement_data$YrResident)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Gender and settlement_data$YrResident
## t = -1.795, df = 409, p-value = 0.07339
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.183570869 0.008389558
## sample estimates:
## cor
## -0.08841148
# Fit regression model
Indonesia_individual_model <- lm(Place_Attachment_Index ~ YrResident +
Age + Gender, data = settlement_data)
summary(Indonesia_individual_model)
##
## Call:
## lm(formula = Place_Attachment_Index ~ YrResident + Age + Gender,
## data = settlement_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.58177 -0.11080 -0.01866 0.08451 1.00455
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.1429070 0.1954449 21.197 <2e-16 ***
## YrResident 0.0020024 0.0013464 1.487 0.138
## Age 0.0004525 0.0028256 0.160 0.873
## Gender -0.2260258 0.1389061 -1.627 0.104
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2094 on 407 degrees of freedom
## Multiple R-squared: 0.01592, Adjusted R-squared: 0.008666
## F-statistic: 2.195 on 3 and 407 DF, p-value: 0.08812
# AIC Value of linear model
AIC(Indonesia_individual_model)
## [1] -112.9781
# Evaluate Model Diagnostics
plot(Indonesia_individual_model)
# Communicate Results
# The results of our multiple linear regression reveal that age, gender, and
# years of residency are not significant predictors of place attachment in
# Indonesia (F(410) = 2.195, adjusted R^2 = 0.009, p = 0.09).
# Null Hypothesis: Place attachment between all MPAs are the same.
# Alternative Hypothesis: Place attachment between at least two
# MPAs are different.
# Examining Data Distributions
ggplot(settlement_data, aes(x = Place_Attachment_Index,
fill = MPAName)) +
geom_histogram() +
facet_grid(.~MPAName) +
labs(x = "Place Attachment Index",
y = "Count",
color = "MPA Name") +
theme_bw() +
theme(legend.position = "bottom",
strip.text.x = element_text(size = 7))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Use Q-Q Plots, skewness, and kurtosis to test for normality
ggplot(settlement_data, aes(sample = Place_Attachment_Index)) +
geom_qq() +
geom_qq_line() +
labs (x = "Theoretical Normal Distribution",
y = "Raw Place Attachment Index Values")
skewness(settlement_data$Place_Attachment_Index)
## [1] 1.254365
kurtosis(settlement_data$Place_Attachment_Index)
## [1] 6.929917
# Non-normal data, so we will use Kruskal-Wallis and Dunn's post-hoc tests
Indonesia_MPA_kruskal_test <- kruskal.test(Place_Attachment_Index ~ MPAName,
data = settlement_data)
Indonesia_MPA_kruskal_test
##
## Kruskal-Wallis rank sum test
##
## data: Place_Attachment_Index by MPAName
## Kruskal-Wallis chi-squared = 16.183, df = 5, p-value = 0.00634
# Results of Kruskal Test are significant, so let's run a Dunn's Test
Indonesia_MPA_dunn_test <- dunn.test(settlement_data$Place_Attachment_Index, settlement_data$MPAName)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 16.1834, df = 5, p-value = 0.01
##
##
## Comparison of x by group
## (No adjustment)
## Col Mean-|
## Row Mean | Kaimana Kofiau d Misool S Selat Da Teluk Ce
## ---------+-------------------------------------------------------
## Kofiau d | -2.688998
## | 0.0036*
## |
## Misool S | -2.118790 1.056939
## | 0.0171* 0.1453
## |
## Selat Da | -0.143330 2.586212 1.987140
## | 0.4430 0.0049* 0.0235*
## |
## Teluk Ce | -2.768019 0.772198 -0.469646 -2.631757
## | 0.0028* 0.2200 0.3193 0.0042*
## |
## Teluk Ma | -2.193081 0.768487 -0.278202 -2.072756 0.113658
## | 0.0142* 0.2211 0.3904 0.0191* 0.4548
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
# Null Hypothesis: Food security between all MPAs are the same.
# Alternative Hypothesis: Food security between at least two
# MPAs are different.
# Examining Data Distributions
ggplot(settlement_data, aes(x = Food_Security_Index,
fill = MPAName)) +
geom_histogram() +
facet_grid(.~MPAName) +
labs(x = "Food Security Index",
y = "Count",
color = "MPA Name") +
theme_bw() +
theme(legend.position = "bottom",
strip.text.x = element_text(size = 7))
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
# Use Q-Q Plots, skewness, and kurtosis to test for normality
ggplot(settlement_data, aes(sample = Food_Security_Index)) +
geom_qq() +
geom_qq_line() +
labs (x = "Theoretical Normal Distribution",
y = "Raw Food Security Index Values")
skewness(settlement_data$Food_Security_Index)
## [1] -0.8122661
kurtosis(settlement_data$Food_Security_Index)
## [1] 4.571809
# Non-normal data, so we will use Kruskal-Wallis and Dunn's post-hoc tests
Indonesia_MPA_kruskal_test2 <- kruskal.test(Food_Security_Index ~ MPAName,
data = settlement_data)
Indonesia_MPA_kruskal_test2
##
## Kruskal-Wallis rank sum test
##
## data: Food_Security_Index by MPAName
## Kruskal-Wallis chi-squared = 13.316, df = 5, p-value = 0.02059
# Results of Kruskal Test are significant, so let's run a Dunn's Test
Indonesia_MPA_dunn_test2 <- dunn.test(settlement_data$Food_Security_Index, settlement_data$MPAName)
## Kruskal-Wallis rank sum test
##
## data: x and group
## Kruskal-Wallis chi-squared = 13.3157, df = 5, p-value = 0.02
##
##
## Comparison of x by group
## (No adjustment)
## Col Mean-|
## Row Mean | Kaimana Kofiau d Misool S Selat Da Teluk Ce
## ---------+-------------------------------------------------------
## Kofiau d | 1.934056
## | 0.0266
## |
## Misool S | 1.843798 -0.503765
## | 0.0326 0.3072
## |
## Selat Da | -0.962296 -2.723294 -2.858532
## | 0.1680 0.0032* 0.0021*
## |
## Teluk Ce | 0.404722 -1.783939 -1.696291 1.505559
## | 0.3428 0.0372 0.0449 0.0661
## |
## Teluk Ma | 1.331516 -0.756317 -0.354337 2.246937 1.114744
## | 0.0915 0.2247 0.3615 0.0123* 0.1325
##
## alpha = 0.05
## Reject Ho if p <= alpha/2
# Make MPA names into a factor with 6 levels
settlement_data2 <- settlement_data %>%
mutate(MPAName = factor(MPAName, levels = c("Kaimana MPA",
"Kofiau dan Pulau Boo MPA",
"Misool Selatan Timur MPA",
"Selat Dampier MPA",
"Teluk Cenderawasih NP",
"Teluk Mayalibit MPA")))
# Define Research Question:
# Do marine tenure, food security index, years of residency, and age
# significantly predict place attachment index in Indonesian households?
# Examine data and possible correlations
# Raw counts
ggplot(data = settlement_data2, aes(x = Household_Marine_Tenure)) +
geom_histogram() +
labs(y = "Count") +
theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = settlement_data2, aes(x = Food_Security_Index)) +
geom_histogram() +
labs(y = "Count") +
theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = settlement_data2, aes(x = YrResident)) +
geom_histogram() +
labs(y = "Count") +
theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = settlement_data2, aes(x = Age)) +
geom_histogram() +
labs(y = "Count") +
theme_bw()
## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.
ggplot(data = settlement_data2, aes(x = Place_Attachment_Index)) +
geom_bar() +
labs(y = "Count") +
theme_bw()
# Possible Correlations
ggplot(data = settlement_data2, aes(x = Household_Marine_Tenure,
y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = Food_Security_Index,
y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = YrResident,
y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = Age, y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
# Multi-collinearity tests
cor.test(settlement_data2$Household_Marine_Tenure,
settlement_data2$Food_Security_Index)
##
## Pearson's product-moment correlation
##
## data: settlement_data2$Household_Marine_Tenure and settlement_data2$Food_Security_Index
## t = -5.8, df = 409, p-value = 1.329e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3627357 -0.1838526
## sample estimates:
## cor
## -0.2756793
cor.test(settlement_data2$Household_Marine_Tenure,
settlement_data2$YrResident)
##
## Pearson's product-moment correlation
##
## data: settlement_data2$Household_Marine_Tenure and settlement_data2$YrResident
## t = 2.036, df = 409, p-value = 0.04239
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.003473309 0.195008737
## sample estimates:
## cor
## 0.1001689
cor.test(settlement_data2$Household_Marine_Tenure,
settlement_data2$Age)
##
## Pearson's product-moment correlation
##
## data: settlement_data2$Household_Marine_Tenure and settlement_data2$Age
## t = -3.0222, df = 409, p-value = 0.002667
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.24107960 -0.05180823
## sample estimates:
## cor
## -0.1477969
cor.test(settlement_data2$Food_Security_Index,
settlement_data2$YrResident)
##
## Pearson's product-moment correlation
##
## data: settlement_data2$Food_Security_Index and settlement_data2$YrResident
## t = 1.5794, df = 409, p-value = 0.115
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.0190111 0.1732856
## sample estimates:
## cor
## 0.07786137
cor.test(settlement_data2$Food_Security_Index,
settlement_data2$Age)
##
## Pearson's product-moment correlation
##
## data: settlement_data2$Food_Security_Index and settlement_data2$Age
## t = 5.4314, df = 409, p-value = 9.616e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1668286 0.3473861
## sample estimates:
## cor
## 0.2593723
cor.test(settlement_data2$YrResident,
settlement_data2$Age)
##
## Pearson's product-moment correlation
##
## data: settlement_data2$YrResident and settlement_data2$Age
## t = 9.553, df = 409, p-value < 2.2e-16
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.3446207 0.5030580
## sample estimates:
## cor
## 0.4271122
# Fit linear regression model without random effect and evaluate model
# diagnostics
place_attachment_lm <- lm(Place_Attachment_Index ~
Household_Marine_Tenure +
Food_Security_Index +
YrResident +
Age, data = settlement_data2)
summary(place_attachment_lm)
##
## Call:
## lm(formula = Place_Attachment_Index ~ Household_Marine_Tenure +
## Food_Security_Index + YrResident + Age, data = settlement_data2)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.57484 -0.11324 -0.00842 0.08892 0.96095
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.719496 0.132809 28.006 < 2e-16 ***
## Household_Marine_Tenure 0.051603 0.017048 3.027 0.00263 **
## Food_Security_Index -0.012266 0.018334 -0.669 0.50383
## YrResident 0.001102 0.001356 0.813 0.41664
## Age 0.004051 0.002806 1.443 0.14969
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2074 on 406 degrees of freedom
## Multiple R-squared: 0.03631, Adjusted R-squared: 0.02682
## F-statistic: 3.825 on 4 and 406 DF, p-value: 0.00459
plot(place_attachment_lm)
# AIC Value of standard linear model
AIC(place_attachment_lm)
## [1] -119.5846
# Fit linear regression model with random effect and evaluate model
# diagnostics
# Re-fit standard linear regression
place_attachment_lm2 <- gls(Place_Attachment_Index ~
Household_Marine_Tenure +
Food_Security_Index +
YrResident +
Age, data = settlement_data2)
# Fit multi-level model with MPA name as random effect
place_attachment_rlm <- lme(Place_Attachment_Index ~
Household_Marine_Tenure +
Food_Security_Index +
YrResident +
Age, random = ~ 1 | MPAName, data = settlement_data2,
na.action=na.exclude)
# Compare AIC Values
AIC(place_attachment_lm2, place_attachment_rlm)
## Warning in AIC.default(place_attachment_lm2, place_attachment_rlm): models are
## not all fitted to the same number of observations
# Plot residuals of multi-level model
plot(place_attachment_rlm)
qqnorm(place_attachment_rlm)
# Model summary
summary(place_attachment_rlm)
## Linear mixed-effects model fit by REML
## Data: settlement_data2
## AIC BIC logLik
## -52.2051 -24.73531 33.10255
##
## Random effects:
## Formula: ~1 | MPAName
## (Intercept) Residual
## StdDev: 0.05852237 0.2094833
##
## Fixed effects: Place_Attachment_Index ~ Household_Marine_Tenure + Food_Security_Index + YrResident + Age
## Value Std.Error DF t-value p-value
## (Intercept) 3.771316 0.15066850 370 25.030557 0.0000
## Household_Marine_Tenure 0.036404 0.02025540 370 1.797269 0.0731
## Food_Security_Index -0.016311 0.01978587 370 -0.824391 0.4102
## YrResident 0.002316 0.00151550 370 1.528322 0.1273
## Age 0.003297 0.00300653 370 1.096582 0.2735
## Correlation:
## (Intr) Hs_M_T Fd_S_I YrRsdn
## Household_Marine_Tenure -0.556
## Food_Security_Index -0.378 0.264
## YrResident 0.139 -0.134 0.015
## Age -0.737 0.174 -0.181 -0.439
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.61399176 -0.63189980 -0.05838216 0.45724255 4.20992865
##
## Number of Observations: 379
## Number of Groups: 5
#Sameer's Logistics Regression Model for determining Place Attachment
#Research Question: Do socioeconomic indicators (Household Material Assets, Household Marine Tenure, Food Security Index and School Enrollment Rate predict binary place attachment to MPAs (Attached/Not Attached)?)
# Possible relationships between covariates and response variable
ggplot(data = settlement_data2, aes(x = Household_Material_Assets,
y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data2, aes(x = Household_Marine_Tenure,
y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = Food_Security_Index,
y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = School_Enrollment_Rate,
y = Place_Attachment_Index)) +
geom_point() +
theme_bw()
#Only School Enrollment Rate shows some relationship with PLace Attachment - slightly higher density of points with a >4.0 place attachment when School Enrollment Rate increases. All other variables show a random distribution of points with respect to Place Attachment
# Look at correlation between covariates and response variable
cor.test(settlement_data$Place_Attachment_Index, settlement_data$Household_Material_Assets)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Place_Attachment_Index and settlement_data$Household_Material_Assets
## t = -1.2861, df = 409, p-value = 0.1991
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.15921755 0.03346906
## sample estimates:
## cor
## -0.0634657
# -0.06
cor.test(settlement_data$Place_Attachment_Index, settlement_data$Household_Marine_Tenure)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Place_Attachment_Index and settlement_data$Household_Marine_Tenure
## t = 3.2469, df = 409, p-value = 0.001263
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.06275387 0.25139540
## sample estimates:
## cor
## 0.1585209
# 0.159
cor.test(settlement_data$Place_Attachment_Index, settlement_data$School_Enrollment_Rate)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Place_Attachment_Index and settlement_data$School_Enrollment_Rate
## t = -2.7473, df = 409, p-value = 0.006275
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.22836305 -0.03837756
## sample estimates:
## cor
## -0.1346072
# -0.135
cor.test(settlement_data$Place_Attachment_Index, settlement_data$Food_Security_Index)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Place_Attachment_Index and settlement_data$Food_Security_Index
## t = -1.0821, df = 409, p-value = 0.2798
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.14938718 0.04352427
## sample estimates:
## cor
## -0.05342996
# -0.053
# Multi-collinearity tests
cor.test(settlement_data$Household_Material_Assets, settlement_data$Household_Marine_Tenure)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Material_Assets and settlement_data$Household_Marine_Tenure
## t = -1.2265, df = 409, p-value = 0.2207
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1563505 0.0364057
## sample estimates:
## cor
## -0.06053677
# correlation of the two variable is statistically non-significant, p-value = 0.221, cor = -0.06
cor.test(settlement_data$Household_Material_Assets, settlement_data$School_Enrollment_Rate)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Material_Assets and settlement_data$School_Enrollment_Rate
## t = 4.1855, df = 409, p-value = 3.486e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1080535 0.2936375
## sample estimates:
## cor
## 0.2026645
# correlation of the two variable is statistically significant, p-value = 3.49e-05, cor = 0.203
cor.test(settlement_data$Household_Material_Assets, settlement_data$Food_Security_Index)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Material_Assets and settlement_data$Food_Security_Index
## t = 3.5123, df = 409, p-value = 0.0004938
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.0756311 0.2634770
## sample estimates:
## cor
## 0.1711086
# correlation of the two variable is statistically significant, p-value = 0.0005, cor = 0.171
cor.test(settlement_data$Household_Marine_Tenure, settlement_data$School_Enrollment_Rate)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Marine_Tenure and settlement_data$School_Enrollment_Rate
## t = -3.3954, df = 409, p-value = 0.0007523
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.25816849 -0.06996526
## sample estimates:
## cor
## -0.165574
# correlation of the two variable is statistically significant, p-value = 0.0008, cor = -0.166
cor.test(settlement_data$Household_Marine_Tenure, settlement_data$Food_Security_Index)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Marine_Tenure and settlement_data$Food_Security_Index
## t = -5.8, df = 409, p-value = 1.329e-08
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.3627357 -0.1838526
## sample estimates:
## cor
## -0.2756793
# correlation of the two variable is statistically significant, p-value = 1.33e-08, cor = -0.276
cor.test(settlement_data$Food_Security_Index, settlement_data$School_Enrollment_Rate)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Food_Security_Index and settlement_data$School_Enrollment_Rate
## t = 2.6431, df = 409, p-value = 0.008531
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.03327849 0.22351807
## sample estimates:
## cor
## 0.1295906
# correlation of the two variable is statistically significant, p-value = 0.009, cor = 0.130
#Statistical significance of correlations between multiple covariate pairs. However, all the correlations are weak (under 0.3) so we proceed with using all the covariates
#Define a function for Rice's Rule
rice_rule <- function(n) {
return(2 * (n^(1/3)))
}
num_bins <- ceiling(rice_rule(length(settlement_data$Place_Attachment_Index)))
#Look at distribution of place attachment values
ggplot(settlement_data, aes(x = Place_Attachment_Index)) +
geom_histogram(bins = num_bins) +
labs(x = "Place Attachment", y = "Frequency", title = "Distribution of Place Attachment values amongst Eastern Indonesian communities") +
theme_bw()
#Distribution is approximately normal with a slight right-skew
mean_PIA <- mean(settlement_data$Place_Attachment_Index)
#As an initial pass, convert Place_Attachment_Index into a binary categorical variable where values below mean are 0 and values above mean are 1. This categorization can be further explored and tweaked
settlement_data <- settlement_data %>%
mutate(PIA_Binary = ifelse(Place_Attachment_Index > mean_PIA, 1, 0))
#Build first logistic regression model with PIA_Binary as response variable
PIA_glm <- glm(PIA_Binary ~ Household_Material_Assets + Household_Marine_Tenure + Food_Security_Index + School_Enrollment_Rate, data = settlement_data)
# Communicate Results
# The results of our logistic regression reveal that Household Material Assets is a significant predictor (p < 0.05) of the binary place attachment variable.
# AIC: 595.94
# household material assets p = 0.035
# household marine tenure p = 0.224
#food security index p = 0.123
# school enrollment rate p = 0.507
#As an initial pass, Household Material Assets is a significant predictor (p < 0.05) of PIA_Binary
#Interesting as Household Material Assets did not show a strong correlation with overall Place Attachment but could be due to binary PIA values.
summary(PIA_glm)
##
## Call:
## glm(formula = PIA_Binary ~ Household_Material_Assets + Household_Marine_Tenure +
## Food_Security_Index + School_Enrollment_Rate, data = settlement_data)
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.46552 0.26301 1.770 0.0775 .
## Household_Material_Assets -0.52772 0.24952 -2.115 0.0350 *
## Household_Marine_Tenure 0.04904 0.04030 1.217 0.2243
## Food_Security_Index 0.06689 0.04324 1.547 0.1227
## School_Enrollment_Rate -0.13420 0.20218 -0.664 0.5072
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## (Dispersion parameter for gaussian family taken to be 0.2454008)
##
## Null deviance: 101.625 on 410 degrees of freedom
## Residual deviance: 99.633 on 406 degrees of freedom
## AIC: 595.94
##
## Number of Fisher Scoring iterations: 2
#No significant outliers from Cook's distance plot
#Normal Q-Q plot shows residuals do not follow a normal distribution with very heavy tails on lower and upper end
#Residuals also do not show a random distribution and homoscedasticity
#Given the AIC value and the residual plots, we will not proceed with this model
plot(PIA_glm)
#Research Question: Do socioeconomic indicators (Household Material Assets, Household Marine Tenure, Food Security Index and School Enrollment Rate predict place attachment to MPAs, when accounting for individual variation of place attachment within MPAs?
# Fit multi-level model with random intercept by MPA as a sensitivity analysis
place_attachment_household_mlm <- lme(Place_Attachment_Index ~
Household_Marine_Tenure +
Household_Material_Assets +
School_Enrollment_Rate +
Food_Security_Index, random = ~ 1 | MPAName, data = settlement_data2, na.action = na.exclude)
summary(place_attachment_household_mlm)
## Linear mixed-effects model fit by REML
## Data: settlement_data2
## AIC BIC logLik
## -68.69349 -41.2237 41.34674
##
## Random effects:
## Formula: ~1 | MPAName
## (Intercept) Residual
## StdDev: 0.0589693 0.2092381
##
## Fixed effects: Place_Attachment_Index ~ Household_Marine_Tenure + Household_Material_Assets + School_Enrollment_Rate + Food_Security_Index
## Value Std.Error DF t-value p-value
## (Intercept) 4.198641 0.13190308 370 31.83125 0.0000
## Household_Marine_Tenure 0.024336 0.02037826 370 1.19422 0.2332
## Household_Material_Assets -0.136370 0.11518354 370 -1.18394 0.2372
## School_Enrollment_Rate -0.194838 0.09097394 370 -2.14169 0.0329
## Food_Security_Index -0.004839 0.01945609 370 -0.24871 0.8037
## Correlation:
## (Intr) Hs_M_T Hs_M_A Sc_E_R
## Household_Marine_Tenure -0.644
## Household_Material_Assets -0.368 0.158
## School_Enrollment_Rate -0.499 0.121 -0.169
## Food_Security_Index -0.553 0.276 -0.041 -0.066
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.82829277 -0.66655643 -0.04793978 0.47614800 4.20203149
##
## Number of Observations: 379
## Number of Groups: 5
# Communicate Results
# The results of our multi-level regression reveal that School Enrollment Rate is a significant predictor (p < 0.05) of place attachment
# AIC: -68.69
#This checks out with our observations of the scatterplots
# household material assets p = 0.237
# household marine tenure p = 0.233
#food security index p = 0.804
# school enrollment rate p = 0.033
#Residuals show mostly homoscedasticity
plot(place_attachment_household_mlm)
#Residuals do not meet normality assumption - heavier tails on both ends
qqnorm(place_attachment_household_mlm)
# Compare AIC Values with logistic regression model - lower AIC value for multi-level model
#Better residual plots and lower AIC value means we choose multi-level model over logistic regression model
AIC(PIA_glm, place_attachment_household_mlm)
## Warning in AIC.default(PIA_glm, place_attachment_household_mlm): models are not
## all fitted to the same number of observations
# Define Research Question:
# Is there a relationship between place attachment index and household material assets,
# household marine tenure, and school enrollment rate?
# Examine data and possible correlations
# Raw Counts
ggplot(data = settlement_data2, aes(x = `Place_Attachment_Index`)) +
geom_histogram(bins = 15) +
theme_bw()
# pretty normally distributed
ggplot(data = settlement_data2, aes(x = `Household_Material_Assets`)) +
geom_histogram(bins = 15) +
theme_bw()
# pretty normally distributed
ggplot(data = settlement_data2, aes(x = `Household_Marine_Tenure`)) +
geom_histogram(bins = 15) +
theme_bw()
# not normally distributed
ggplot(data = settlement_data, aes(x = `School_Enrollment_Rate`)) +
geom_histogram(bins = 15) +
theme_bw()
# pretty skewed
# Possible Correlations
ggplot(data = settlement_data, aes(x = `Household_Material_Assets`, y = `Place_Attachment_Index`)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = `Household_Marine_Tenure`, y = `Place_Attachment_Index`)) +
geom_point() +
theme_bw()
ggplot(data = settlement_data, aes(x = `School_Enrollment_Rate`, y = `Place_Attachment_Index`)) +
geom_point() +
theme_bw()
#unclear if any of those have any sort of correlation
# Numerically view correlation
cor.test(settlement_data$Place_Attachment_Index, settlement_data$`Household_Material_Assets`)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Place_Attachment_Index and settlement_data$Household_Material_Assets
## t = -1.2861, df = 409, p-value = 0.1991
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.15921755 0.03346906
## sample estimates:
## cor
## -0.0634657
# very low.... -0.06
cor.test(settlement_data$Place_Attachment_Index, settlement_data$`Household_Marine_Tenure`)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Place_Attachment_Index and settlement_data$Household_Marine_Tenure
## t = 3.2469, df = 409, p-value = 0.001263
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.06275387 0.25139540
## sample estimates:
## cor
## 0.1585209
# slightly better .. 0.158
cor.test(settlement_data$Place_Attachment_Index, settlement_data$`School_Enrollment_Rate`)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Place_Attachment_Index and settlement_data$School_Enrollment_Rate
## t = -2.7473, df = 409, p-value = 0.006275
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.22836305 -0.03837756
## sample estimates:
## cor
## -0.1346072
# still low ... -0.135
# Multi-collinearity tests
cor.test(settlement_data$`Household_Material_Assets`, settlement_data$`Household_Marine_Tenure`)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Material_Assets and settlement_data$Household_Marine_Tenure
## t = -1.2265, df = 409, p-value = 0.2207
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.1563505 0.0364057
## sample estimates:
## cor
## -0.06053677
# pass
cor.test(settlement_data$`Household_Material_Assets`, settlement_data$`School_Enrollment_Rate`)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Material_Assets and settlement_data$School_Enrollment_Rate
## t = 4.1855, df = 409, p-value = 3.486e-05
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## 0.1080535 0.2936375
## sample estimates:
## cor
## 0.2026645
# pass
cor.test(settlement_data$`Household_Marine_Tenure`, settlement_data$`School_Enrollment_Rate`)
##
## Pearson's product-moment correlation
##
## data: settlement_data$Household_Marine_Tenure and settlement_data$School_Enrollment_Rate
## t = -3.3954, df = 409, p-value = 0.0007523
## alternative hypothesis: true correlation is not equal to 0
## 95 percent confidence interval:
## -0.25816849 -0.06996526
## sample estimates:
## cor
## -0.165574
# pass
# ok so no strong correlations, not worried about co-linearity either
# try to fit model and see what happens
# Fit regression model
Indonesia_household_model <- lm(Place_Attachment_Index ~ `Household_Material_Assets` +
`Household_Marine_Tenure` + `School_Enrollment_Rate`,
data = settlement_data)
summary(Indonesia_household_model)
##
## Call:
## lm(formula = Place_Attachment_Index ~ Household_Material_Assets +
## Household_Marine_Tenure + School_Enrollment_Rate, data = settlement_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63066 -0.12002 -0.00832 0.09008 0.94731
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.08103 0.09169 44.509 < 2e-16 ***
## Household_Material_Assets -0.07023 0.10312 -0.681 0.49622
## Household_Marine_Tenure 0.04589 0.01626 2.822 0.00501 **
## School_Enrollment_Rate -0.17568 0.08433 -2.083 0.03784 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.207 on 407 degrees of freedom
## Multiple R-squared: 0.0383, Adjusted R-squared: 0.03121
## F-statistic: 5.403 on 3 and 407 DF, p-value: 0.001182
# Evaluate Model Diagnostics
plot(Indonesia_household_model)
# plots look fine, no outliers
# Communicate Results
# The results of our multiple linear regression reveal that household marine tenure,
# and school enrollment rate are significant predictors of place attachment in
# Indonesia (F(407) = 5.403, adjusted R^2 = 0.03121, p = 0.0012).
# household material assets t = 0.496
# household marine tenure t = 0.005 **
# school enrollment rate t = 0.0378 *
# This model is weak, let's see what it looks like without household material assets
# Fit regression model
Indonesia_household_model2 <- lm(Place_Attachment_Index ~ `Household_Marine_Tenure`
+ `School_Enrollment_Rate`,
data = settlement_data)
summary(Indonesia_household_model2)
##
## Call:
## lm(formula = Place_Attachment_Index ~ Household_Marine_Tenure +
## School_Enrollment_Rate, data = settlement_data)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.63033 -0.12315 -0.00653 0.09096 0.93604
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 4.05556 0.08366 48.477 < 2e-16 ***
## Household_Marine_Tenure 0.04620 0.01625 2.844 0.00468 **
## School_Enrollment_Rate -0.18692 0.08264 -2.262 0.02423 *
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.2068 on 408 degrees of freedom
## Multiple R-squared: 0.0372, Adjusted R-squared: 0.03248
## F-statistic: 7.882 on 2 and 408 DF, p-value: 0.0004377
# AIC Value
AIC(Indonesia_household_model2)
## [1] -123.964
# Evaluate Model Diagnostics
plot(Indonesia_household_model2)
# plots look fine
# The results of our multiple linear regression reveal that household marine tenure,
# and school enrollment rate are significant predictors of place attachment in
# Indonesia (F(408) = 7.882, adjusted R^2 = 0.03248, p = 0.0004).
# household marine tenure t = 0.00468 **
# school enrollment rate t = 0.0242 *
# better, but not sure what this tells us. maybe we try a multilevel model.
# Trying a mixed effects model with different predictor variables. A mix of individual and
# household variables to try and build a more significant model.
# Try a multilevel/mixed effects model using vars we know have a correlation
# with place attachment
# Fit multi-level model with random intercept by plot
place_attachment_mem <- lme(Place_Attachment_Index ~
School_Enrollment_Rate +
Household_Marine_Tenure +
YrResident,
random = ~ 1 | MPAName, data = settlement_data2,
na.action=na.exclude)
# Compare AIC Values between all mixed effect models so far
AIC(place_attachment_mem, place_attachment_rlm, place_attachment_household_mlm)
## Warning in AIC.default(place_attachment_mem, place_attachment_rlm,
## place_attachment_household_mlm): models are not all fitted to the same number
## of observations
# mlm has the lowest AIC, with the mem having second lowest
# Plot residuals of multi-level model
plot(place_attachment_mem)
qqnorm(place_attachment_mem)
# Model summary
summary(place_attachment_mem)
## Linear mixed-effects model fit by REML
## Data: settlement_data2
## AIC BIC logLik
## -72.75938 -49.19782 42.37969
##
## Random effects:
## Formula: ~1 | MPAName
## (Intercept) Residual
## StdDev: 0.06266477 0.2074717
##
## Fixed effects: Place_Attachment_Index ~ School_Enrollment_Rate + Household_Marine_Tenure + YrResident
## Value Std.Error DF t-value p-value
## (Intercept) 4.052218 0.10060027 371 40.28039 0.0000
## School_Enrollment_Rate -0.243205 0.08936277 371 -2.72154 0.0068
## Household_Marine_Tenure 0.025942 0.01921014 371 1.35042 0.1777
## YrResident 0.003459 0.00135611 371 2.55101 0.0111
## Correlation:
## (Intr) Sc_E_R Hs_M_T
## School_Enrollment_Rate -0.760
## Household_Marine_Tenure -0.567 0.188
## YrResident -0.268 -0.121 -0.067
##
## Standardized Within-Group Residuals:
## Min Q1 Med Q3 Max
## -2.85107406 -0.64953050 -0.03015959 0.48822683 4.14707607
##
## Number of Observations: 379
## Number of Groups: 5
# Make table with AIC results
# Create a data frame with model names and AIC values
aic_model1 <- AIC(place_attachment_rlm)
aic_model2 <- AIC(place_attachment_household_mlm)
aic_model3 <- AIC(place_attachment_mem)
aic_table <- data.frame(
Model = c("Mixed Effect Model to Determine Place Attachment Attempt 1",
"Mixed Effect Model to Determine Place Attachment Attempt 2",
"Mixed Effect Model to Determine Place Attachment Attempt 3"),
AIC = c(aic_model1,
aic_model2,
aic_model3)
)
# Print the table
print(aic_table)
## Model AIC
## 1 Mixed Effect Model to Determine Place Attachment Attempt 1 -52.20510
## 2 Mixed Effect Model to Determine Place Attachment Attempt 2 -68.69349
## 3 Mixed Effect Model to Determine Place Attachment Attempt 3 -72.75938
# make pretty using gt
aic_table_gt <- aic_table %>%
gt() %>%
tab_header(
title = "AIC Comparison of Mixed-Effects Models"
) %>%
fmt_number(
columns = vars(AIC),
decimals = 2
) %>%
cols_label(
Model = "Model",
AIC = "AIC Value"
) %>%
tab_spanner_delim(
delim = " ",
columns = everything()
)
## Warning: Since gt v0.3.0, `columns = vars(...)` has been deprecated.
## • Please use `columns = c(...)` instead.
# Print the table
aic_table_gt
| AIC Comparison of Mixed-Effects Models | |
| Model | AIC Value |
|---|---|
| Mixed Effect Model to Determine Place Attachment Attempt 1 | −52.21 |
| Mixed Effect Model to Determine Place Attachment Attempt 2 | −68.69 |
| Mixed Effect Model to Determine Place Attachment Attempt 3 | −72.76 |
# save the table
gtsave(aic_table_gt, "aic_table.png")
## file:////var/folders/1k/wfyvwjys5cs9vp3xc7rpf18h0000gn/T//RtmptnazJi/file1ba124b510d4.html screenshot completed
# Summary Table
indonesia_summary <- settlement_data %>%
group_by(MPAName) %>%
summarise(count = n(),
mean_place_attachment = mean(Place_Attachment_Index,
na.rm = TRUE),
median_place_attachment = median(Place_Attachment_Index, na.rm = TRUE),
variance_place_attachment = var(Place_Attachment_Index, na.rm = TRUE),
sd_place_attachment = sd(Place_Attachment_Index, na.rm = TRUE)) %>%
ungroup()
# Rename columns
colnames(indonesia_summary) <- c("MPA Name", "Number of Observations",
"Mean Place Attachment",
"Median Place Attachment",
"Variance Place Attachment",
"SD Place Attachment")
# Export Table
write.csv(indonesia_summary, "Indonesia_Summary_Statistics.csv")
indonesia_gt <- indonesia_summary %>% gt() %>%
fmt_number(columns = c(`Mean Place Attachment`,
`Median Place Attachment`,
`Variance Place Attachment`,
`SD Place Attachment`),
decimals = 2)
indonesia_gt
| MPA Name | Number of Observations | Mean Place Attachment | Median Place Attachment | Variance Place Attachment | SD Place Attachment |
|---|---|---|---|---|---|
| Kaimana MPA | 65 | 3.96 | 3.96 | 0.03 | 0.17 |
| Kofiau dan Pulau Boo MPA | 32 | 4.04 | 4.03 | 0.01 | 0.10 |
| Misool Selatan Timur MPA | 76 | 3.99 | 4.02 | 0.02 | 0.16 |
| Selat Dampier MPA | 67 | 3.97 | 3.95 | 0.03 | 0.18 |
| Teluk Cenderawasih NP | 119 | 4.02 | 4.02 | 0.03 | 0.18 |
| Teluk Mayalibit MPA | 52 | 4.14 | 4.00 | 0.14 | 0.37 |
gtsave(indonesia_gt,
"indonesia_summary_table.png",
path = here())
## file:////var/folders/1k/wfyvwjys5cs9vp3xc7rpf18h0000gn/T//RtmptnazJi/file1ba127cd7c74.html screenshot completed
# Create mean of other variables table
indonesia_mean_summary <- settlement_data %>%
summarise(mean_household_material_assets = mean(Household_Material_Assets,
na.rm = TRUE),
mean_household_marine_tenure = mean(Household_Marine_Tenure,
na.rm = TRUE),
mean_food_security_index = mean(Food_Security_Index,
na.rm = TRUE),
mean_place_attachment_index = mean(Place_Attachment_Index,
na.rm = TRUE),
mean_school_enrollment_rate = mean(School_Enrollment_Rate,
na.rm = TRUE),
mean_years_residency = mean(YrResident, na.rm = TRUE),
mean_household_age = mean(Age, na.rm = TRUE),
mean_gender = mean(Gender, na.rm = TRUE))
# Rename columns
colnames(indonesia_mean_summary) <- c("Household Material Assets",
"Household Marine Tenure",
"Food Security Index",
"Place Attachment Index",
"School Enrollment Rate",
"Years of Residency",
"Age",
"Gender")
indonesia_mean_summary <- indonesia_mean_summary %>%
pivot_longer(cols = `Household Material Assets`:Gender,
names_to = "Variable Name",
values_to = "Mean")
# Export Table
indonesia_gt2 <- indonesia_mean_summary %>% gt() %>%
fmt_number(columns = c(Mean),
decimals = 2)
indonesia_gt2
| Variable Name | Mean |
|---|---|
| Household Material Assets | 0.48 |
| Household Marine Tenure | 2.40 |
| Food Security Index | 3.41 |
| Place Attachment Index | 4.01 |
| School Enrollment Rate | 0.81 |
| Years of Residency | 29.06 |
| Age | 44.69 |
| Gender | 0.91 |
gtsave(indonesia_gt2,
"indonesia_mean_summary_table.png",
path = here())
## file:////var/folders/1k/wfyvwjys5cs9vp3xc7rpf18h0000gn/T//RtmptnazJi/file1ba123a5b816.html screenshot completed